      program main
      implicit none
      integer*4,parameter:: Nl=5,Nwmin=-100,Nwmax=100,Nk=256
      double precision pi,dw,tn,tnn,tnnn,dk,eps,zrn,xmu,Akw,Akw1,Akw2,
     &                 qx,qy,
     &                 disp(Nl,-Nk:2*Nk),dispQ(Nl,-Nk:2*Nk),
     &                 dSC(Nl,-Nk:2*Nk),AFM(Nl),
     &                 Vhyb(Nl-1),D0(Nl),E0(Nl),Hamil0(Nl,Nl),
     &                 Amax,Hamil(4*Nl,4*Nl),
     &                 Eval(4*Nl),RWORK(16*Nl)
      complex*16 xi,com,cdum(Nl,Nl),cdum2(4*Nl,4*Nl),
     &              WORK(Nl*Nl),WORK2(16*Nl*Nl)
      integer*4 i,j,k,kx,ky,ii,jj,kk,kmax,kf,kshift,kinit,kend,iflag,
     &          INFO,IPIV(Nl),IPIV2(4*Nl)
      character*13 filename
      
      pi=acos(-1.0d0)
      xi=cmplx(0.0d0,1.0d0)

      open(unit=1,file='input.txt',status='old')
      read(1,*)tn    !energy scale
      read(1,*)tnn,tnnn,xmu
      read(1,*)(Vhyb(j),j=1,Nl-1)  !Between nearest layers
      read(1,*)(E0(j),j=1,Nl)
      read(1,*)(D0(j),j=1,Nl)
      read(1,*)(AFM(j),j=1,Nl)
      read(1,*)eps,dw
      read(1,*)kshift
      read(1,*)zrn
      close(1)
      tnn=tnn*tn
      tnnn=tnnn*tn
      xmu=xmu*tn
      Vhyb(:)=Vhyb(:)*tn
      D0(:)=D0(:)*tn
      E0(:)=E0(:)*tn
      AFM(:)=AFM(:)*tn

      dk=2.0d0*pi/(2*Nk)

      do i=1,Nl
      do k=0,Nk
        kx=k
        ky=k-kshift
        disp(i,k)=-2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &            -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &            -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
        dispQ(i,k)=2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &            -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &            -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
        dSC(i,k)=0.5d0*D0(i)*(cos(kx*dk)-cos(ky*dk))
      end do
      end do

c      Akw=0.0d0
c      do k=Nk,kshift,-1   ! determine kf for inner FS
c        do i=1,Nl
c          Hamil0(i,i)=disp(i,k)  !underlying FS
c        end do
c        do i=1,Nl-1
c          Hamil0(i,i+1)=Vhyb(i)
c          Hamil0(i+1,i)=Vhyb(i)
c        end do
c
c        com=xi*eps   !omega=0
c
c        cdum(:,:)=-Hamil0(:,:)
c        do i=1,Nl
c          cdum(i,i)=cdum(i,i)+com
c        end do
c
c        call ZGETRF(Nl,Nl,cdum,Nl,IPIV,INFO)
c        call ZGETRI(Nl,cdum,Nl,IPIV,WORK,Nl*Nl,INFO)
c
c        Akw1=Akw
c
c        Akw=0.0d0
c        do i=1,Nl
c          Akw=Akw-aimag(cdum(i,i))
c        end do
c        Akw=Akw/pi
c
c        if(Akw.lt.Akw1)then
c          kf=k+1
c          exit
c        end if
c      end do

c      Amax=0.0d0 
c      do k=kf,kshift,-1   ! determine kf for outer FS
c        do i=1,Nl
c          Hamil0(i,i)=disp(i,k)  !underlying FS
c        end do
c        do i=1,Nl-1
c          Hamil0(i,i+1)=Vhyb(i)
c          Hamil0(i+1,i)=Vhyb(i)
c        end do
c
c        com=xi*eps   !omega=0
c
c        cdum(:,:)=-Hamil0(:,:)
c        do i=1,Nl
c          cdum(i,i)=cdum(i,i)+com
c        end do
c
c        call ZGETRF(Nl,Nl,cdum,Nl,IPIV,INFO)
c        call ZGETRI(Nl,cdum,Nl,IPIV,WORK,Nl*Nl,INFO)
c
c        Akw=0.0d0
c        do i=1,Nl
c          Akw=Akw-aimag(cdum(i,i))
c        end do
c        Akw=Akw/pi
c
c        if(Akw.gt.Amax)then
c          kmax=k
c          Amax=Akw
c        end if 
c      end do

      kinit=Nk/2-Nk/4
      kend=Nk/2+Nk/4

      open(unit=10,file='Akw_00-PP.txt',status='unknown')
      open(unit=20,file='Evec_00-PP.txt',status='unknown')
      do k=kinit,kend
        kk=k
        if(kk.gt.Nk)kk=2*Nk-k

        Hamil(:,:)=0.0d0
        do i=1,Nl
          Hamil(i,i)=       disp(i,kk)
          Hamil(i+Nl,i+Nl)=-disp(i,kk)
          Hamil(i,i+Nl)=dSC(i,kk)
          Hamil(i+Nl,i)=dSC(i,kk)

          Hamil(i+2*Nl,i+2*Nl)= dispQ(i,kk)
          Hamil(i+3*Nl,i+3*Nl)=-dispQ(i,kk)
          Hamil(i+2*Nl,i+3*Nl)=-dSC(i,kk)
          Hamil(i+3*Nl,i+2*Nl)=-dSC(i,kk)
        end do
        do i=1,Nl-1
          Hamil(i,i+1)=Vhyb(i)
          Hamil(i+1,i)=Vhyb(i)
          Hamil(i+Nl,i+1+Nl)=-Vhyb(i)
          Hamil(i+1+Nl,i+Nl)=-Vhyb(i)

          Hamil(i+2*Nl,i+1+2*Nl)=Vhyb(i)
          Hamil(i+1+2*Nl,i+2*Nl)=Vhyb(i)
          Hamil(i+3*Nl,i+1+3*Nl)=-Vhyb(i)
          Hamil(i+1+3*Nl,i+3*Nl)=-Vhyb(i)
        end do

        do i=1,Nl
          Hamil(i,i+2*Nl)=AFM(i)
          Hamil(i+2*Nl,i)=AFM(i)
          Hamil(i+Nl,i+3*Nl)=AFM(i)
          Hamil(i+3*Nl,i+Nl)=AFM(i)
        end do

        do ii=Nwmin,Nwmax
          com=ii*dw+xi*eps

          cdum2(:,:)=-Hamil(:,:)
          do i=1,4*Nl
            cdum2(i,i)=cdum2(i,i)+com
          end do

          call ZGETRF(4*Nl,4*Nl,cdum2,4*Nl,IPIV2,INFO)
          call ZGETRI(4*Nl,cdum2,4*Nl,IPIV2,WORK2,16*Nl*Nl,INFO)

          Akw=0.0d0
          do i=1,Nl
            Akw=Akw-aimag(cdum2(i,i))
          end do
          Akw=Akw/pi

          write(10,'(3F14.6)')k*dk*sqrt(2.0d0)/pi,ii*dw,Akw

          if(ii.eq.0)write(100,*)k*dk/pi,Akw
          if(ii.eq.0)Akw1=Akw
        end do
        write(10,*)
        if(Akw1.gt.20.0d0)then
          call DSYEV('V','U',4*Nl,Hamil,4*Nl,Eval,RWORK,16*Nl,INFO)
          do ii=1,4*Nl
            if(abs(Eval(ii)).lt.0.01d0)then
              write(20,'(3F14.6)') k*dk*sqrt(2.0d0)/pi,Akw1,Eval(ii)
              write(20,'(5F14.6)')(Hamil(j     ,ii),j=1,Nl)
              write(20,'(5F14.6)')(Hamil(j+  Nl,ii),j=1,Nl)
              write(20,'(5F14.6)')(Hamil(j+2*Nl,ii),j=1,Nl)
              write(20,'(5F14.6)')(Hamil(j+3*Nl,ii),j=1,Nl)
              write(20,*)
            end if
          end do
        end if
      end do
      close(10)
      close(20)

c----------------------------------
      do jj=0,8
      kshift=Nk/16*jj
      kinit=Nk/2-Nk/4+kshift
      kend=Nk/2+Nk/4+kshift


      write(filename,'(a8,I1,a4)')'Akw_diag',jj,'.txt'
      open(unit=10,file=filename,status='unknown')
      do k=kinit,kend
        kx=k
        ky=k-kshift*2
        do i=1,Nl
          disp(i,k)=-2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &              -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &              -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
          dispQ(i,k)=2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &              -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &              -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
          dSC(i,k)=0.5d0*D0(i)*(cos(kx*dk)-cos(ky*dk))
        end do

        do i=1,Nl
          Hamil(i,i)=       disp(i,k)
          Hamil(i+Nl,i+Nl)=-disp(i,k)
          Hamil(i,i+Nl)=dSC(i,k)
          Hamil(i+Nl,i)=dSC(i,k)

          Hamil(i+2*Nl,i+2*Nl)= dispQ(i,k)
          Hamil(i+3*Nl,i+3*Nl)=-dispQ(i,k)
          Hamil(i+2*Nl,i+3*Nl)=-dSC(i,k)
          Hamil(i+3*Nl,i+2*Nl)=-dSC(i,k)
        end do
        do i=1,Nl-1
          Hamil(i,i+1)=Vhyb(i)
          Hamil(i+1,i)=Vhyb(i)
          Hamil(i+Nl,i+1+Nl)=-Vhyb(i)
          Hamil(i+1+Nl,i+Nl)=-Vhyb(i)

          Hamil(i+2*Nl,i+1+2*Nl)=Vhyb(i)
          Hamil(i+1+2*Nl,i+2*Nl)=Vhyb(i)
          Hamil(i+3*Nl,i+1+3*Nl)=-Vhyb(i)
          Hamil(i+1+3*Nl,i+3*Nl)=-Vhyb(i)
        end do

        do i=1,Nl
          Hamil(i,i+2*Nl)=AFM(i)
          Hamil(i+2*Nl,i)=AFM(i)
          Hamil(i+Nl,i+3*Nl)=AFM(i)
          Hamil(i+3*Nl,i+Nl)=AFM(i)
        end do

        do ii=Nwmin,Nwmax
          com=ii*dw+xi*eps

          cdum2(:,:)=-Hamil(:,:)
          do i=1,4*Nl
            cdum2(i,i)=cdum2(i,i)+com
          end do

          call ZGETRF(4*Nl,4*Nl,cdum2,4*Nl,IPIV2,INFO)
          call ZGETRI(4*Nl,cdum2,4*Nl,IPIV2,WORK2,16*Nl*Nl,INFO)

          Akw=0.0d0
          do i=1,Nl
            Akw=Akw-aimag(cdum2(i,i))
          end do
          Akw=Akw/pi

          write(10,'(4F14.6)')kx*dk,ky*dk,ii*dw,Akw

        end do
        write(10,*)
      end do
      close(10)

      end do

c----------------------------------
      do jj=0,8
      kshift=Nk/2+Nk/16*jj
      kinit=-Nk
      kend=Nk


      write(filename,'(a8,I1,a4)')'Akw_xfix',jj,'.txt'
      open(unit=10,file=filename,status='unknown')
      do k=kinit,kend
        kx=kshift
        ky=k
        do i=1,Nl
          disp(i,k)=-2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &              -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &              -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
          dispQ(i,k)=2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &              -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &              -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
          dSC(i,k)=0.5d0*D0(i)*(cos(kx*dk)-cos(ky*dk))
        end do

        do i=1,Nl
          Hamil(i,i)=       disp(i,k)
          Hamil(i+Nl,i+Nl)=-disp(i,k)
          Hamil(i,i+Nl)=dSC(i,k)
          Hamil(i+Nl,i)=dSC(i,k)

          Hamil(i+2*Nl,i+2*Nl)= dispQ(i,k)
          Hamil(i+3*Nl,i+3*Nl)=-dispQ(i,k)
          Hamil(i+2*Nl,i+3*Nl)=-dSC(i,k)
          Hamil(i+3*Nl,i+2*Nl)=-dSC(i,k)
        end do
        do i=1,Nl-1
          Hamil(i,i+1)=Vhyb(i)
          Hamil(i+1,i)=Vhyb(i)
          Hamil(i+Nl,i+1+Nl)=-Vhyb(i)
          Hamil(i+1+Nl,i+Nl)=-Vhyb(i)

          Hamil(i+2*Nl,i+1+2*Nl)=Vhyb(i)
          Hamil(i+1+2*Nl,i+2*Nl)=Vhyb(i)
          Hamil(i+3*Nl,i+1+3*Nl)=-Vhyb(i)
          Hamil(i+1+3*Nl,i+3*Nl)=-Vhyb(i)
        end do

        do i=1,Nl
          Hamil(i,i+2*Nl)=AFM(i)
          Hamil(i+2*Nl,i)=AFM(i)
          Hamil(i+Nl,i+3*Nl)=AFM(i)
          Hamil(i+3*Nl,i+Nl)=AFM(i)
        end do

        do ii=Nwmin,Nwmax
          com=ii*dw+xi*eps

          cdum2(:,:)=-Hamil(:,:)
          do i=1,4*Nl
            cdum2(i,i)=cdum2(i,i)+com
          end do

          call ZGETRF(4*Nl,4*Nl,cdum2,4*Nl,IPIV2,INFO)
          call ZGETRI(4*Nl,cdum2,4*Nl,IPIV2,WORK2,16*Nl*Nl,INFO)

          Akw=0.0d0
          do i=1,Nl
            Akw=Akw-aimag(cdum2(i,i))
          end do
          Akw=Akw/pi

          write(10,'(4F14.6)')kx*dk,ky*dk,ii*dw,Akw

        end do
        write(10,*)
      end do
      close(10)

      end do
c        ----------------------------------------
      do i=1,Nl
      do k=0,Nk
        kx=Nk-k
        ky=k
        disp(i,k)=-2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &            -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &            -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
        dispQ(i,k)=2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &            -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(i)-xmu
     &            -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
        dSC(i,k)=0.5d0*D0(i)*(cos(kx*dk)-cos(ky*dk))
      end do
      end do
      kinit=0
      kend=Nk/2
      Amax=0.0d0

      iflag=0
      open(unit=10,file='Akw_P0-0P.txt',status='unknown')
      do k=kinit,kend
        kk=k
        if(kk.gt.Nk)kk=2*Nk-k

        do i=1,Nl
          Hamil(i,i)=       disp(i,kk)
          Hamil(i+Nl,i+Nl)=-disp(i,kk)
          Hamil(i,i+Nl)=dSC(i,kk)
          Hamil(i+Nl,i)=dSC(i,kk)

          Hamil(i+2*Nl,i+2*Nl)= dispQ(i,kk)
          Hamil(i+3*Nl,i+3*Nl)=-dispQ(i,kk)
          Hamil(i+2*Nl,i+3*Nl)=-dSC(i,kk)
          Hamil(i+3*Nl,i+2*Nl)=-dSC(i,kk)
        end do
        do i=1,Nl-1
          Hamil(i,i+1)=Vhyb(i)
          Hamil(i+1,i)=Vhyb(i)
          Hamil(i+Nl,i+1+Nl)=-Vhyb(i)
          Hamil(i+1+Nl,i+Nl)=-Vhyb(i)

          Hamil(i+2*Nl,i+1+2*Nl)=Vhyb(i)
          Hamil(i+1+2*Nl,i+2*Nl)=Vhyb(i)
          Hamil(i+3*Nl,i+1+3*Nl)=-Vhyb(i)
          Hamil(i+1+3*Nl,i+3*Nl)=-Vhyb(i)
        end do

        do i=1,Nl
          Hamil(i,i+2*Nl)=AFM(i)
          Hamil(i+2*Nl,i)=AFM(i)
          Hamil(i+Nl,i+3*Nl)=AFM(i)
          Hamil(i+3*Nl,i+Nl)=AFM(i)
        end do

        do ii=Nwmin,Nwmax
          com=ii*dw+xi*eps

          cdum2(:,:)=-Hamil(:,:)
          do i=1,4*Nl
            cdum2(i,i)=cdum2(i,i)+com
          end do

          call ZGETRF(4*Nl,4*Nl,cdum2,4*Nl,IPIV2,INFO)
          call ZGETRI(4*Nl,cdum2,4*Nl,IPIV2,WORK2,16*Nl*Nl,INFO)

          Akw=0.0d0
          do i=1,Nl
            Akw=Akw-aimag(cdum2(i,i))
          end do
          Akw=Akw/pi
          write(10,'(3F14.6)')k*dk*sqrt(2.0d0)/pi,ii*dw,Akw

          if(ii.eq.0)then
            if(kk.gt.Nk/4)then
              if(Akw.gt.Amax.and.iflag.eq.0)then
                kmax=kk
                Amax=Akw
                !print*,kmax,Amax
              else
                iflag=1
              end if
            end if 
          end if
        end do
        write(10,*)
      end do
      close(10)

      open(unit=12,file='Akw_kfinner.txt',status='unknown')
        qy=kmax*dk
        qx=pi-qy
        disp(:,0)=-2.0d0*tn*(cos(qx)+cos(qy))
     &        -4.0d0*tnn*cos(qx)*cos(qy)+E0(:)-xmu
     &        -2.0d0*tnnn*(cos(2*qx)+cos(2*qy))
        dispQ(:,0)=2.0d0*tn*(cos(qx)+cos(qy))
     &        -4.0d0*tnn*cos(qx)*cos(qy)+E0(:)-xmu
     &        -2.0d0*tnnn*(cos(2*qx)+cos(2*qy))
        dSC(:,0)=0.5d0*D0(:)*(cos(qx)-cos(qy))

        do i=1,Nl
          Hamil(i,i)=       disp(i,0)
          Hamil(i+Nl,i+Nl)=-disp(i,0)
          Hamil(i,i+Nl)=dSC(i,0)
          Hamil(i+Nl,i)=dSC(i,0)

          Hamil(i+2*Nl,i+2*Nl)= dispQ(i,0)
          Hamil(i+3*Nl,i+3*Nl)=-dispQ(i,0)
          Hamil(i+2*Nl,i+3*Nl)=-dSC(i,0)
          Hamil(i+3*Nl,i+2*Nl)=-dSC(i,0)
        end do
        do i=1,Nl-1
          Hamil(i,i+1)=Vhyb(i)
          Hamil(i+1,i)=Vhyb(i)
          Hamil(i+Nl,i+1+Nl)=-Vhyb(i)
          Hamil(i+1+Nl,i+Nl)=-Vhyb(i)

          Hamil(i+2*Nl,i+1+2*Nl)=Vhyb(i)
          Hamil(i+1+2*Nl,i+2*Nl)=Vhyb(i)
          Hamil(i+3*Nl,i+1+3*Nl)=-Vhyb(i)
          Hamil(i+1+3*Nl,i+3*Nl)=-Vhyb(i)
        end do

        do i=1,Nl
          Hamil(i,i+2*Nl)=AFM(i)
          Hamil(i+2*Nl,i)=AFM(i)
          Hamil(i+Nl,i+3*Nl)=AFM(i)
          Hamil(i+3*Nl,i+Nl)=AFM(i)
        end do

        write(12,*)'# kf_inner=',qx,qy
        do ii=Nwmin,Nwmax
          com=ii*dw+xi*eps

          cdum2(:,:)=-Hamil(:,:)
          do i=1,4*Nl
            cdum2(i,i)=cdum2(i,i)+com
          end do

          call ZGETRF(4*Nl,4*Nl,cdum2,4*Nl,IPIV2,INFO)
          call ZGETRI(4*Nl,cdum2,4*Nl,IPIV2,WORK2,16*Nl*Nl,INFO)

          Akw=0.0d0
          do i=1,Nl
            Akw=Akw-aimag(cdum2(i,i))
          end do
          Akw=Akw/pi

          write(12,'(2F14.6)')ii*dw,Akw
        end do
      close(12)

      open(unit=11,file='Akw0.txt',status='unknown')
      do kx=0,Nk
      do ky=0,Nk
        disp(:,0)=-2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &        -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(:)-xmu
     &        -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
        dispQ(:,0)=2.0d0*tn*(cos(kx*dk)+cos(ky*dk))
     &        -4.0d0*tnn*cos(kx*dk)*cos(ky*dk)+E0(:)-xmu
     &        -2.0d0*tnnn*(cos(2*kx*dk)+cos(2*ky*dk))
        dSC(:,0)=0.5d0*D0(:)*(cos(kx*dk)-cos(ky*dk))

        do i=1,Nl
          Hamil(i,i)=       disp(i,0)
          Hamil(i+Nl,i+Nl)=-disp(i,0)
          Hamil(i,i+Nl)=dSC(i,0)
          Hamil(i+Nl,i)=dSC(i,0)

          Hamil(i+2*Nl,i+2*Nl)= dispQ(i,0)
          Hamil(i+3*Nl,i+3*Nl)=-dispQ(i,0)
          Hamil(i+2*Nl,i+3*Nl)=-dSC(i,0)
          Hamil(i+3*Nl,i+2*Nl)=-dSC(i,0)
        end do
        do i=1,Nl-1
          Hamil(i,i+1)=Vhyb(i)
          Hamil(i+1,i)=Vhyb(i)
          Hamil(i+Nl,i+1+Nl)=-Vhyb(i)
          Hamil(i+1+Nl,i+Nl)=-Vhyb(i)

          Hamil(i+2*Nl,i+1+2*Nl)=Vhyb(i)
          Hamil(i+1+2*Nl,i+2*Nl)=Vhyb(i)
          Hamil(i+3*Nl,i+1+3*Nl)=-Vhyb(i)
          Hamil(i+1+3*Nl,i+3*Nl)=-Vhyb(i)
        end do

        do i=1,Nl
          Hamil(i,i+2*Nl)=AFM(i)
          Hamil(i+2*Nl,i)=AFM(i)
          Hamil(i+Nl,i+3*Nl)=AFM(i)
          Hamil(i+3*Nl,i+Nl)=AFM(i)
        end do

        cdum2(:,:)=-Hamil(:,:)
        do i=1,4*Nl
          cdum2(i,i)=cdum2(i,i)+xi*eps
        end do

        call ZGETRF(4*Nl,4*Nl,cdum2,4*Nl,IPIV2,INFO)
        call ZGETRI(4*Nl,cdum2,4*Nl,IPIV2,WORK2,16*Nl*Nl,INFO)

        Akw=0.0d0
        do i=1,Nl
          Akw=Akw-aimag(cdum2(i,i))
        end do
        Akw=Akw/pi

        write(11,'(3F14.6)')kx*dk,ky*dk,Akw
      end do
      write(11,*)
      end do
      end
